function [ blip, tout ] = get_blip_math_lsim(amiga_model, ... 
    audio_filter_mode, filter_method, length)

if (nargin < 3)
	filter_method = 'lsim';
end

if (nargin < 4)
    length = 'short';
end

filter_short_len = get_filter_short_len();
fs = get_sampling_fq();
spectrum_cut = get_spectrum_cut();

[ts1, ts2] = get_amiga_filter_ts(amiga_model, audio_filter_mode);

if strcmpi(length, 'full') || strcmpi(filter_method, 'fq')
    tout = linspace(0, (fs-2)/fs, fs-1);
elseif strcmpi(length, 'short')
    tout = linspace(0, (filter_short_len-1) / fs, filter_short_len);
end

if strcmpi(filter_method, 'lsim')
    blip = firwin(amiga_model, fs, spectrum_cut, tout, false);
else
    blip = tout*0;
    blip(1) = 1;
end

if strcmpi(filter_method, 'lsim')
    blip = fir_min_phase(blip);
end

blip = lsim(ts1*ts2, blip', tout')';

if ~strcmpi(filter_method, 'lsim')
    if ~strcmpi(filter_method, 'none')
        blip = fir_min_phase(blip);
    end
end

blip = filter_impulse(amiga_model, blip, filter_method);

if strcmpi(length, 'short') && strcmpi(filter_method, 'fq')
    tout = tout(1:1:filter_short_len);
    blip = blip(1:1:filter_short_len);
end

blip = blip / trapz(blip);

end